Refinement  of  Seafloor  Elevation 
Using  Acoustic  Backscatter 

Andrew  E.  Johnson  and  Martial  Hebert 


CMU-RI  TR-95-01 


Robotics  Institute 
Carnegie  Mellon  University 


©  1995  Carnegie  Mellon  University 


This  research  was  supported  by  the  National  Science  Foundation  under  contract  DCR- 
8604199. 

The  views  and  conclusions  contained  in  this  document  are  those  of  the  authors  and  should 
not  be  interpreted  as  representing  the  official  policies,  either  expressed  or  implied,  of  the 
U.S.  government. 


19950425  042 


Abstract 


We  propose  an  algorithm  for  the  reconstruction  of  elevation  and  material  prop¬ 
erty  maps  of  the  seafloor  using  a  side  scan  sonar  backscatter  image  and  sparse 
bathymetric  points  co-registered  within  the  image.  Given  a  path  for  the  sensor, 
the  reconstruction  is  corrected  for  the  movement  of  the  fish  during  the  image 
generation  process.  To  perform  reconstruction,  an  arbitrary  but  computable 
scattering  model  is  assumed  for  the  seafloor  backscatter. 

The  algorithm  uses  the  sparse  bathymetric  data  to  generate  an  initial  estimate 
for  the  elevation  map  which  is  then  iteratively  refined  to  fit  the  backscatter 
image  by  minimizing  a  global  error  functional.  Concurrently,  the  parameters 
of  the  scattering  model  are  determined  on  a  coarse  grid  in  the  image  by  fitting 
the  assumed  scattering  model  to  the  backscatter  data.  The  elevation  surface 
and  the  scattering  parameter  maps  converge  to  their  best  fit  shape  and  values 
given  the  backscatter  data.  The  reconstruction  is  corrected  for  the  movement  of 
the  sensor  by  initially  doing  local  reconstructions  in  sensor  coordinates  and 
then  transforming  the  local  reconstructions  to  a  global  coordinate  system  and 
performing  the  reconstruction  again. 

The  algorithm  supports  different  scattering  models,  so  it  can  be  applied  to  dif¬ 
ferent  underwater  environments  and  sonar  sensors.  In  addition  to  the  elevation 
map  of  the  seafloor,  the  parameters  of  the  scattering  model  at  every  point  in  the 
image  are  generated.  Since  these  parameters  describe  material  properties  of 
the  seafloor,  the  maps  of  the  scattering  model  parameters  can  be  used  to  seg¬ 
ment  the  seafloor  by  material  type. 
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1  Introduction 

The  undersea  world  is  an  unfriendly  environment  for  navigation  sensors  mainly  because  the 
sensors  must  sense  through  the  water  medium.  Water  severely  attenuates  and  distorts  light,  so 
optical  sensors  produce  poor  results  at  far  range.  Sonar,  on  the  other  hand,  is  more  suitable  for 
sensing  underwater.  Water  does  not  attenuate  sound  energy  nearly  as  much  as  light  energy,  so 
sonar  has  a  greater  range.  Sonar  is  also  less  susceptible  to  refraction.  These  benefits  come  at 
the  price  of  lower  resolution  and  increased  noisiness  of  the  sonar  data.  However,  given  the 
task  of  mapping  and  navigation  underwater,  sonar  is  the  best  alternative. 

Oeanographers  and  geologists  have  a  strong  desire  to  map  the  seafloor  because  the  terrain  and 
material  composition  of  the  seafloor  hold  clues  to  understanding  the  geological  processes  that 
transform  the  earth.  Scattering  parameters  are  related  to  actual  physical  properties  of  the 
ocean  floor,  so  measuring  and  mapping  them  will  help  geologists  understand  the  physical 
makeup  of  the  ocean  floor  thus  establishing  useful  patterns  that  can  be  related  to  the  evolution 
of  the  planet.  Furthermore,  these  physical  properties  also  give  clues  to  underlying  materials  in 
the  seafloor  which  can  be  useful  for  mining  and  rare  element  acquisition.  Elevation  maps  have 
obvious  geological  significance,  and  they  can  be  used  for  autonomous  underwater  vehicle 
(AUV)  navigation  because  an  AUV  must  know  the  shape  of  the  seafloor  around  it  in  order  to 
navigate  close  to  the  ocean  bottom. 


1.1  Sidescan  Sonar 

A  simple  sonar  sensor  consists  of  one  cylindrical  source  that  creates  a  conic  acoustic  beam 
pattern  that  is  symmetric  around  the  axis  of  the  source.  This  type  of  sonar  will  measure  the 
range  to  the  first  surface  it  encounters  within  the  cone  and  the  intensity  or  echo  of  the  return, 
however,  the  position  of  the  surface  cannot  be  localized  within  the  cone.  The  angle  of  the  cone 
is  to  large  to  make  range  sensing  with  these  type  of  simple  sensors  feasible.  Side  scan  sonars 
are  arrays  of  sonars  arrange  in  a  straight  line.  When  the  sensor  generates  an  acoustic  pulse 
(pings),  the  pulses  from  the  individual  sonars  are  phased  in  such  a  way  that  acoustic  beam  pat¬ 
terns  from  each  individual  sonar  interfere  to  shrink  the  overall  beam  pattern  of  the  array  in  the 
direction  perpendicular  to  the  line  of  the  sensors.  Simply  put,  the  array  acts  like  an  acoustic 
diffraction  grating.  When  the  array  is  pinged,  an  acoustic  pulse  is  generated  that  travels 
through  the  water.  When  the  pulse  encounters  the  ocean  floor  part  of  the  pulse  is  scattered 
back  to  the  sensor.  The  sensor  measures  these  echoes  from  the  seafloor  for  a  short  time  inter¬ 
val  after  the  ping.  Hence,  given  that  the  speed  of  sound  in  water  is  known,  recording  the  mea¬ 
surements  of  the  sensor  will  generate  acoustic  return  as  a  function  of  range  from  the  sensor. 

Side  scan  sonar  are  the  sensors  that  are  generally  used  for  mapping  of  the  seafloor.  The  sonar 
array  is  mounted  on  a  torpedo  shaped  platform  (sonar  fish)  that  is  towed  through  the  water  by 
a  surface  ship.  An  acoustic  image  is  generated  by  pulling  the  platform  through  the  water;  the 
acoustic  return  vs.  range  function  for  each  ping  constitute  the  rows  of  the  image  which  are 
stacked  to  form  the  image.  If  the  fish  is  towed  in  a  horizontal  line,  then  the  image  is  acoustic 
return  as  a  function  of  range  from  the  sensor  and  position  of  the  sensor  along  the  line  S(R,y). 


1 


FIGURE  1.  A  typical  unfiltered  backscatter  image  of  the  seafloor 


1.2  Definition  of  Acoustic  Backscatter 

Acoustic  Backscatter  is  the  ratio  of  sound  energy  impinging  on  a  surface  to  that  leaving  the 
surface.  Generally,  it  is  a  function  of  incidence  angle  of  the  acoustic  pulse  and  material  prop¬ 
erties  of  the  surface.  Because  backscatter  contains  geometric  information  about  the  surface 
being  imaged,  it  is  the  quantity  desired  when  trying  to  reconstruct  the  seafloor  elevation  map. 
Unfortunately,  the  acoustic  backscatter  from  the  surface  is  not  directly  measured  by  the  sen¬ 
sor.  Many  other  processes  relating  to  the  transmission  of  sound  through  the  water  and  the 
equipment  used  to  measure  it  affect  the  data  collection.  The  sonar  equation  presented  in 
Equation  1  details  these  processes  [33]. 

lOlogS  =  RL-SL-aR- TVG  +  2  TL  -  lOlogA  -  RBP  -  SBP  (EQ  1) 

The  backscatter  S  from  the  seafloor  is  a  function  of  receiver  level  RL,  source  level  SL,  water 
column  attenuation  which  is  a  function  of  range  aR,  the  time  varying  gain  of  the  sensor  TVG, 
the  transmission  loss  of  the  acoustic  pulse  due  to  spherical  spreading  TL,  the  equivalent 
ensonified  area  A,  receiver  beam  pattern  RBP  and  the  source  beam  pattern  SBP.  If  all  quanti¬ 
ties  except  for  receiver  level  can  be  independently  measured,  then  backscatter  can  be  obtained 
from  the  receiver  level  measured  by  the  sensor.  An  image  of  acoustic  backscatter  from  a  rela¬ 
tively  flat  seabed  is  given  in  Figure  1.  Sidescan  sonar  calibration,  which  requires  complex 
testing  facilities  and  detailed  modeling  of  sonar  sensors,  has  been  studied  extensively 
[20]  [23]  [30]  [31]. 

1.3  Problem  Statement 

Side  scan  sonar  data  in  the  form  of  a  backscatter  image  does  not  directly  convey  the  elevation 
of  the  seafloor.  However,  given  knowledge  of  the  physics  of  sonar  scattering  from  the  water- 
seafloor  interface  (i.e.,  a  known  scattering  model)  it  is  possible,  if  appropriate  assumptions  are 
made,  to  convert  this  backscatter  image  into  an  elevation  map  of  the  ocean  bottom.  The  scat¬ 
tering  model  with  respect  to  sound  of  the  seafloor  along  with  its  correct  parameters  are 
inverted  to  obtain  the  incidence  angle  of  the  sonar  pulse.  These  incidence  angles  are  then  con¬ 
verted  to  elevation  derivatives  which  are  integrated  together  to  form  the  best  estimate  of  the 
surface  of  the  seafloor  given  the  backscatter  data.  The  problem  our  algorithm  seeks  to  address 
is  the  following: 

Given  an  arbitrary  but  computable  backscatter  model,  generate  accurate  and 
detailed  global  elevation  and  scattering  model  parameter  maps  of  the  seafloor 
using  a  sidescan  sonar  backscatter  image,  sparse  (co-registered)  bathymetric 
data  and  the  path  of  the  sensor  during  the  image  generation  process. 
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There  has  been  a  great  deal  of  research  conducted  to  determine  a  general  and  physically 
meaningful  scattering  model  for  the  seafloor  [8][13][28].  To  date  the  most  detailed  model  for 
high  frequency  acoustic  backscatter  combines  models  for  three  different  scattering  processes 
from  the  seafloor  (volume  sediment,  Kirchoff  approximation  and  composite  roughness)  into 
one  scattering  model  [16]  [24].  However,  the  most  detailed  models  do  not  assume  that  mea¬ 
surements  of  local  surface  shape  of  the  seafloor  are  available,  so  the  models  are  presented  as 
backscatter  as  a  function  of  grazing  angle,  not  incidence  angle.  However,  there  exist  some 
recent  work  that  combines  known  bathymetry  to  generate  local  slope  estimates  and  backscat¬ 
ter  to  estimate  scattering  model  parameters  of  the  seafloor  [11][21][25].  These  approaches  use 
dense  bathymetric  images,  and,  as  a  result,  are  restricted  to  data  gathered  from  bathymetric 
sidescan  sonars. 

Existing  methods  of  seafloor  reconstruction  have  used  a  specific  scattering  model  and  some¬ 
what  arbitrarily  selected  scattering  model  parameters  for  the  seafloor  to  determine  incidence 
angles  [19].  If  the  model  and  parameters  are  not  guessed  correctly,  then  the  integration  will 
produce  a  surface  that  is  not  correct.  Furthermore,  if  the  integration  does  not  have  fixed  start¬ 
ing  and  ending  points,  accumulation  of  errors  will  force  the  surface  away  from  its  real  value. 
Our  previous  work  [17]  did  not  use  arbitrarily  set  scattering  parameters,  but  it  did  not  account 
for  the  deviation  of  the  path  of  the  fish  from  a  straight  line;  hence  it  could  not  generate  accu¬ 
rate  global  elevation  reconstructions. 


2  Detailed  Problem  Statement 

Ideally  we  would  like  to  generate  elevation  maps  of  the  seafloor  using  only  sonar  intensity 
images.  However,  this  is  not  possible  because  determining  elevation  from  intensity  is  a  highly 
underconstrained  problem.  Echo  intensity  is  related  to  the  angle  of  incidence  a  of  the  sonar 
pulse  with  the  surface  through  the  acoustic  scattering  model  a  at  the  interface  between  the 
water  and  the  seafloor.  This  scattering  model  is  a  function  of  many  scattering  model  parame¬ 
ters  p  and  the  angle  of  incidence.  Furthermore,  the  angle  of  incidence  of  the  surface  is  a  func¬ 
tion  of  the  elevation  z,  its  derivatives  zx  and  zy  and  the  three  dimensional  position  of  the 
sensor  at  each  point.  If  we  parameterize  the  surface  in  terms  of  its  horizontal  coordinates  (x,y) 
we  have  the  following  relationship  for  all  intensity  values  S: 

s  =  o(a(zx,  zy,  z,  (X,  Y,  Z,  0,  <p,  co)  sensor)  ,p)  (EQ  2) 

Thus,  the  scattering  model  has  ten  or  more  variables  depending  on  the  number  of  scattering 
model  parameters.  Determining  the  elevation  at  each  point  can  be  achieved  by  solving 
Equation  2,  for  z,  but  this  cannot  be  done  without  knowledge  of  all  the  other  variables.  To 
make  matters  worse,  we  do  not  know  the  functional  form  of  the  scattering  model  with  respect 
to  its  parameters.  Assumptions  must  be  made  and  additional  data  must  be  added  to  invert  this 
equation 

To  solve  this  problem,  we  have  added  two  additional  sets  of  data.  First,  we  assume  that  we 
know  the  position  of  the  sensor  for  all  intensity  values  collected.  Second,  we  assume  that  we 
have  some  sparse  elevation  data  points  co-registered  within  the  sonar  intensity  image.  The 
position  of  the  fish  can  be  obtained  from  position  sensors  on  the  sonar  fish  (pressure  depth, 
altimeter,  transponder  nets,  compass,  etc.).  The  co-registered  elevation  data  points  can  be 
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obtained  from  feature  matching  between  the  intensity  image  and  some  known  global  elevation 
map,  acoustic  beacons,  bathymetric  side  scan  sonars,  etc. 

We  will  assume  that,  for  each  data  set  being  used,  we  know  the  functional  form  of  the  scatter¬ 
ing  model  appropriate  for  the  seafloor  being  imaged  and  the  sonar  sensor  being  used.  The 
algorithm  does  not  depend  on  a  particular  form  for  the  scattering  model;  the  only  requirement 
is  that  the  model  and  its  derivatives  with  respect  to  its  parameters  have  a  computable  form  for 
all  angles  of  incidence. 

Given  these  additions  and  assumptions,  we  now  know  the  exact  form  of  the  scattering  model, 
where  the  data  was  collected,  and  we  have  a  vague  idea  of  the  shape  of  the  seafloor.  However, 
we  do  not  know  any  of  the  scattering  model  parameters  and  we  do  not  know  the  elevation  and 
its  derivatives  at  each  point;  these  are  the  unknowns  that  our  algorithm  will  determine. 

The  basic  idea  behind  the  algorithm  will  be  to  use  the  sparse  elevation  points  (hereafter  called 
control  points)  to  generate  an  initial  estimate  for  the  elevation  map  and  its  derivatives.  Then 
these  initial  estimates  will  be  iteratively  adjusted  to  create  a  surface  that  fits  the  intensity  data. 
At  each  iteration  the  parameters  of  the  scattering  model  will  be  estimated  based  on  the  current 
shape  of  the  seafloor.  The  end  results  are  elevation  and  scattering  parameter  maps  which  can 
be  used  for  mapping  and  navigation. 

The  path  of  the  fish  as  it  is  towed  through  the  water  affects  the  acoustic  returns  from  the  seaf¬ 
loor  because  the  backscatter  from  the  seafloor  is  a  function  of  fish  pose.  Furthermore,  devia¬ 
tions  of  the  path  from  a  straight  line  can  cause  the  sensor  to  scan  the  same  region  of  the 
seafloor  from  two  different  pings,  thus  causing  overlap  in  the  backscatter  image.  However,  if 
the  path  of  the  fish  is  known,  the  pose  of  the  fish  for  each  ping  can  be  used  to  correct  for  the 
movement  of  the  fish  during  the  acoustic  image  generation.  As  a  result,  the  reconstruction  will 
have  a  local  phase  where  the  fish  is  assumed  to  move  in  a  straight  line  and  the  cylindrical 
geometry  of  the  image  is  maintained,  and  a  global  phase  that  corrects  for  the  movement  of  the 
fish  by  transforming  the  intensity  data  and  the  local  elevation  estimates  to  cartesian  global 
coordinates. 

The  problem  then  will  be  to  determine  the  elevation  map  and  scattering  parameter  maps  of  the 
seafloor  given  a  sidescan  sonar  image,  sparse  elevation  control  points,  the  pose  of  the  fish  for 
all  of  the  data  collected  and  the  functional  form  of  the  scattering  model  used. 

3  Algorithmic  Overview 

Before  going  over  the  algorithm  it  is  useful  to  define  some  coordinate  systems  used  in  the 
reconstruction. 

Fish  coordinates:  The  coordinate  system  attached  to  the  fish  for  the  current 
ping.  This  coordinate  system  is  cylindrical  in  that  its  coordinates  are  (R,y,z). 

This  coordinate  system  is  different  for  each  ping  because  the  position  of  the 
fish  changes  for  each  ping. 

Global  coordinates:  The  cartesian  coordinate  system  to  which  all  of  the  local 
reconstructions  are  transformed.  It  is  fixed  to  the  position  of  the  fish  at  the  first 
ping  in  the  image. 
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FIGURE  2.  (A)  Geometry  of  a  side  scan  sonar  and  (B)  image  generation  coordinate  systems 


A  visualization  of  these  coordinate  systems  is  given  in  Figure  2.  Assuming  that  locally  the  fish 
moves  in  a  straight  horizontal  line,  the  elevation  along  each  ping  can  be  determined  in  fish 
coordinates  for  that  ping  independently  of  the  position  of  the  fish.  Given  this  assumption,  the 
reconstruction  can  proceed  in  two  stages.  In  the  first  stage  the  surface  is  reconstructed  locally 
for  each  ping  only  using  data  from  nearby  pings.  This  will  give  a  correct  reconstruction  of  the 
surface  in  fish  coordinates  for  that  ping.  The  second  stage  converts  all  the  reconstructions  to 
global  coordinates,  merges  overlapping  data  and  completes  by  performing  a  global  recon¬ 
struction  in  cartesian  coordinates.  The  following  steps  detail  the  stages  in  the  algorithm. 

1  Generate  Initial  Estimates 

To  start  the  algorithm,  an  initial  estimate  of  the  elevation  surface  must  be  created  from  the 
sparse  elevation  control  points.  This  initial  estimate  will  ensure  that  the  starting  point  for  the 
surface  is  close  to  its  actual  value.  If  the  assumption  is  made  that  locally  the  fish  moves  per¬ 
pendicularly  to  the  swath  and  in  a  horizontal  straight  line,  then  locally  the  derivatives  of  the 
elevation  surface  can  be  calculated  by  taking  finite  differences  between  consecutive  rows  and 
columns.  The  end  result  of  this  first  interpolation  are  three  maps  embedded  in  the  echo  inten¬ 
sity  image  that  correspond  to  the  elevation  at  each  pixel  z(p,s),  the  derivative  of  the  surface 
with  respect  to  the  range  zs(p,s)  and  the  derivative  of  the  surface  with  respect  to  distance  trav¬ 
elled  zp(p,s). 

2  Create  Local  Maps 

Once  the  initial  estimates  for  the  local  derivatives  have  been  determined  local  reconstructions 
can  commence.  For  each  ping,  subimages  of  intensity,  elevation  and  its  derivatives  are  made 
from  the  parts  of  the  images  and  maps  that  surround  the  current  ping.  These  are  the  local  maps 
for  the  current  ping. 

3  Estimate  Local  Scattering  Parameters  with  Non-linear  Minimization 

The  next  step  in  the  local  reconstruction  is  to  determine  the  scattering  model  parameters  at 
each  point  in  the  image.  To  do  this,  we  make  the  justifiable  assumption  that  pixels  that  are 
close  to  each  other  in  the  image  have  similar  scattering  parameters.  Thus,  given  this  grouping 
plus  the  fact  that  we  know  the  elevation  and  its  derivatives  and  hence  that  angle  of  incidence 
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FIGURE  3.  Block  diagram  for  map  local  generation 


at  each  point  (locally  we  are  ignoring  fish  pose),  we  can  do  a  non-  linear  fit  of  the  scattering 
model  to  the  data  (S,a)  to  determine  the  scattering  parameters. 

4  Reconstruct  Elevation  Map 

Once  the  estimates  of  the  scattering  parameters  have  been  determined,  we  are  in  a  position  to 
refine  the  initial  estimate  of  the  elevation  and  its  derivatives  given  the  dense  intensity  data. 
The  idea  behind  this  stage  is  to  create  a  global  error  functional  that  measures  the  divergence  of 
the  surface  from  the  data.  Minimizing  this  functional  will  force  the  surface  into  a  shape  that  is 
consistent  with  the  intensity  data.  Procedures  of  this  sort  are  very  common  in  computer  vision 
and  are  commonly  called  shape  from  shading  techniques  [14],  Numerical  stability  and  integra- 
bility  of  the  surface  require  the  addition  of  smoothness  terms  to  the  functional  that  enforce 
integrability  of  the  surface. 

5  Iterate  on  Local  Maps 

One  simple  pass  through  these  procedures  is  not  sufficient  to  determine  the  surface  or  the  scat¬ 
tering  parameters  because  they  are  inheritly  coupled.  However,  iterating  by  repeatedly  esti¬ 
mating  the  parameters  and  then  refining  the  surface  will  cause  the  surface  to  converge. 

This  procedure  will  generate  local  elevation  and  scattering  parameter  maps.  Creating  these 
local  maps  for  each  ping  will  generate  a  good  estimate  of  the  elevation  of  the  seafloor  with 
respect  to  the  local  fish  coordinates  for  that  ping. 

Figure  3  shows  the  block  diagram  for  the  local  algorithm.  First,  sparse  bathymetric  control 
points  and  a  given  scattering  model  are  used  to  estimate  the  scattering  model  parameters  and 
elevation  map  every  where  in  the  image.  Then,  dense  intensity  data  is  used  to  iteratively  refine 
the  estimates  of  the  scattering  parameters  and  the  elevation  map.  Iteration  continues  until  the 
elevation  values  and  scattering  parameters  have  converged  every  where  in  the  image. 

6  Transform  Local  Reconstructions  to  Global  Coordinates 

The  local  reconstructions  are  transformed  to  global  coordinates  using  the  position  of  the  fish. 
It  is  common  for  multiple  data  points  to  fall  into  the  same  grid  cell  in  the  global  map  so  spe¬ 
cial  attention  must  be  paid  to  combining  these  points  so  that  a  minimal  amount  of  information 
is  lost.  The  general  idea  is  to  average  the  elevation  values  while  keeping  all  of  the  intensity 
data  around  so  that  it  can  be  used  in  the  reconstruction. 


7  Repeat  Local  Algorithm  on  Global  Maps 

The  global  algorithm  is  essentially  the  same  as  the  local  algorithm.  First  initial  estimates  of 
the  seafloor  elevation  and  its  derivatives  with  respect  to  x  and  y  are  determined.  Then  the  scat¬ 
tering  parameters  are  estimated  in  regions  using  nonlinear  fitting  of  the  scattering  model  to  all 
of  the  intensity  data.  Then  the  elevation  surface  is  refined  using  all  of  the  scattering  data  and 
the  initial  estimates  for  the  elevation.  The  latter  two  stages  are  repeated  until  the  surface 
changes  very  little. 


4  Sensor  and  Data  Specifications 

We  used  real  and  synthetic  data  to  test  and  debug  our  reconstruction  algorithm.  The  real  data 
was  collected  by  members  of  the  Deep  Submergence  Laboratory  at  Woods  Hole  Oceano¬ 
graphic  Institution  in  1991  during  there  multisensor  survey  of  the  Juan  de  Fuca  Ridge.  The 
data  sets  we  received  came  from  two  different  geological  provinces  (a  ridge  flank  and  an  axial 
valley).  Simple  sets  of  synthetic  data  were  generated  to  test  our  algorithm  for  strengths  and 
weaknesses. 

4.1  Sensor 

The  particular  sensor  from  which  the  data  was  acquired  was  a  120-kHz  phase-difference  side- 
scan  sonar  developed  jointly  by  the  Deep  Submergence  Laboratory  at  Woods  Hole  Oceano¬ 
graphic  Institution,  the  Applied  Physics  Laboratory  at  the  University  of  Washington,  and 
Acoustic  Marine  Systems,  Inc.  of  Redmond  Washington  [31].  There  are  two  sidescan  sonar 
arrays  on  each  side  of  the  sonar  fish.  This  sensor  is  able  to  generate  acoustic  images  by  aver¬ 
aging  the  acoustic  receiver  levels  from  both  of  the  arrays.  It  is  also  able  to  generate  dense  co¬ 
registered  bathymetry  of  the  seafloor  by  taking  the  phase  difference  of  the  acoustic  signals 
between  the  two  arrays.  The  phase  difference  can  be  converted  directly  into  seafloor  elevation 
in  fish  coordinates  given  the  phase  difference  and  the  range  to  the  surface.  The  geometry  asso¬ 
ciated  with  this  sensor  is  given  in  Figure  2.  The  incidence  angle  a,  which  is  a  key  factor  in 
determining  the  acoustic  backscatter  from  the  surface,  is  the  angle  between  the  vector  from 
the  sensor  to  the  surface  R  and  the  surface  normal  n.  The  grazing  angle  0  is  the  angle  between 
R  and  the  horizontal  and  is  calculated  from  the  phase  between  the  two  arrays.  Calibrated  mea¬ 
surements  of  source  level,  transmission  beam  pattern,  receiver  beam  pattern,  system  gains  and 
transmission  losses  through  the  medium  have  been  made  for  this  sensor,  so  the  receiver  level 
at  the  sensor  can  be  converted  directly  into  acoustic  backscatter  from  the  seafloor.  Also  avail¬ 
able  are  full  six  dimensional  pose  measurements  for  the  fish  which  come  from  positional  sen¬ 
sors  on  the  fish  and  an  acoustic  transponder  net. 

This  sensor  provides  the  necessary  information  to  test  our  algorithm  on  real  data.  The  six 
dimensional  positional  information  makes  the  transformation  from  fish  coordinates  to  world 
coordinates  possible.  The  calibration  of  the  sensor  makes  available  backscatter  images  of  the 
seafloor.  The  bathymetric  map  provides  the  elevation  control  points  as  well  as  a  ground  truth 
for  verifying  the  reconstructed  elevation  map.  For  the  purposes  of  testing  the  algorithm  on  a 
sparse  set  of  elevation  values,  only  discrete  elevation  control  points  were  used  from  the  bathy¬ 
metric  image.  It  should  be  noted  that  sensors  that  generate  dense  co-registered  bathymetric 
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Filtered  and  Unfiltered  Backscatter 


Filtered  and  Unfiltered  Bathymetry 


sample 
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FIGURE  4.  Acoustic  backscatter  and  bathymetry  for  a  typical  sonar  ping  before  and  after  filtering 


maps  of  the  seafloor  are  not  common,  so  algorithms  for  determining  seafloor  elevation  from 
acoustic  images  are  still  extremely  useful. 


4.2  Filtering 

The  images  generated  using  the  AMS/DSL  120  sensor  are  extremely  noisy  due  to  multiple 
echoes,  the  statistical  nature  of  acoustic  backscatter  and  noise  in  the  sensors,  so  filtering  of  the 
data  was  required.  First,  a  3x3  median  filter  was  applied  to  the  backscatter  and  bathymetric 
images  to  eliminate  isolated  spurious  measurements.  Then,  a  modified  median  filter  with  vari¬ 
able  window  size  that  eliminates  only  dropouts  in  the  data  was  applied.  The  thought  behind 
this  filtering  procedure  is  that  dropouts  are  much  more  likely  to  be  due  to  noise  while  positive 
spikes  are  likely  to  correspond  to  interesting  features  in  the  image.  The  final  filter  applied  was 
a  graduated  non  convexity  filter  which  smooths  data  while  it  preserves  discontinuities.  The 
combination  of  these  filters  provided  smooth  data  yet  did  not  eliminate  any  of  the  interesting 
features  in  the  image.  Before  and  after  images  of  acoustic  backscatter  and  bathymetry  are 
given  in  Figure  4. 

The  pose  of  the  fish  also  required  filtering  because  of  noise  in  the  sensors  and  low  resolution 
of  some  of  the  measurements  (especially  global  X  and  Y).  In  general,  a  median  filter  of  vari¬ 
able  window  size  and  then  a  gaussian  filter  were  applied  to  all  of  the  measurements  indepen- 


dently.  Filtering  of  the  pose  considering  the  dependence  of  the  parameters  was  attempted,  but 
was  later  deemed  unnecessary. 


4.3  Applicable  Scattering  Models 

The  algorithm  was  tested  using  many  different  scattering  models  which  are  listed  in  Table  1. 
The  models  range  from  simple,  as  in  the  case  of  the  Lambertian  model  which  has  one  param¬ 
eter,  to  complex,  as  in  the  case  of  the  composite  roughness  model  which  has  four  parameters. 
The  complex  models  give  a  more  accurate  description  of  back  scatter  from  the  seafloor  at  the 
expense  of  being  difficult  to  compute  and  use.  Even  more  complex  models  can  be  made  by 
adding  models  together  or  blending  them  so  that  they  are  applied  only  at  the  range  of  angles 
for  which  they  are  physically  feasible  [24].  The  lambertian  scattering  model  is  the  simplest 
and  most  widely  used  scattering  model  in  acoustics  and  computer  vision.  The  Torrance  and 
Sparrow  model  and  the  diffuse  model  are  optical  reflection  models  commonly  used  in  com¬ 
puter  vision.  In  these  three  scattering  models,  each  parameter  has  a  physical  meaning:  the 
albedo  p  determines  the  ability  of  the  seafloor  to  absorb  sound,  C  and  e  are  measurements  of 
the  roughness  of  the  seafloor  surface  and  a  is  a  constant  that  determines  the  specularity  of  the 
seafloor.  The  Kirchoff  approximation,  composite  roughness  and  sediment  volume  models  are 
common  scattering  models  from  underwater  acoustics.  The  Kirchoff  approximation  is  a  well 
known  model  that  is  applicable  near  normal  incidence  (a  =  0°)  of  the  acoustic  pulse  on  the 
surface.  The  sediment  volume  scattering  model  gives  the  functional  form  of  the  acoustic 
backscatter  that  originates  from  the  acoustic  pulse  penetrating  the  seafloor  surface  and  then 
being  scattered  back  by  a  perfectly  flat  homogeneous  surface.  The  composite  roughness 
model  explains  the  backscatter  caused  by  the  small  scale  roughness  of  the  surface  and  it  is 
mainly  applicable  at  moderate  incidence  angles.  The  parameters  of  these  models,  which  are 
explained  in  Table  1,  relate  to  the  different  physical  properties  of  the  materials  being  imaged 
by  the  sonar  sensor. 

Figure  5  shows  the  scattering  models  with  the  best  fits  to  the  data  collected  from  an  underwa¬ 
ter  ridge  flank.  The  best  parameters  of  the  models  were  determined  through  the  use  of  a 
genetic  algorithm  for  determining  functional  parameters  of  complex  and  possibly  non-differ- 
entiable  functions.  The  scattering  model  with  the  least  fit  error  is  the  complex  model  pre¬ 
sented  by  Mourad  and  Jackson  [24]  which  is  a  continuous  combination  of  the  composite 
roughness,  volume,  and  Kirchoff  approximations  models. 

4.4  Synthetic  Data 

In  addition  to  the  real  data,  sets  of  synthetic  data  were  generated  to  test  the  algorithm  in  a  con¬ 
trolled  environment.  This  required  simulation  of  the  fish  traveling  on  a  well  known  path  while 
scanning  an  object  whose  shape,  scattering  model  and  scattering  model  parameters  were  well 
known.  Using  data  sets  like  this  allowed  us  to  test  the  efficacy  of  our  algorithm  in  estimating 
the  shape  of  the  surface,  estimating  the  scattering  parameters  given  the  correct  scattering 
model,  and  merging  multiple  views  of  the  same  object.  For  most  of  the  examples  using  syn¬ 
thetic  data  in  this  paper,  the  fish  moves  in  a  semi-circular  path  while  scanning  a  cylindrical 
surface  with  a  lambertian  scattering  model.  This  path  and  shape  of  the  object  are  simple 
enough  to  make  analysis  of  the  results  of  the  algorithm  straightforward,  yet  at  the  same  time 
complex  enough  to  fully  test  all  components  of  the  algorithm.  Gaussian  noise  whose  magni- 
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TABLE  1.  Scattering  models  and  description  of  their  parameters  that  were  tested 


Name 

Form 

Parameters 

Lambertian 

o  =  pcosa 

p  =  albedo 

Diffuse3 

a  = 

f  .  r 

3^cosa  +  esina 

p  =  albedo 
e  =  roughness 

Torrance  & 
Sparrowb 

/ 

a  =  pcosa 

V 

(1  -C)  +2  Cexp 

-a 

„  2 
,2a  JJ 

p  =  albedo 

a  =  rms  surface  roughness 

C  =  weight  of  body  vs.  surface 
scattering 

Kirchoff 

Approximation3 

a  =  — 

871 

bqcR902 

f  .  4r  4V 

^sina  +  cosa  1 

_ 

2  r 

+  r 

27 

-2  r 

r  =  y/2-1 ;  y  =  surface  roughness 
spectrum  exponent 

P  =  surface  roughness  spectral 
strength 

p  =ratio  of  sediment  mass  den¬ 
sity  to  water  mass  density 

v  =  ratio  of  sediment  sound 
speed  to  water  sound  speed 

ka  is  the  acoustic  pulse  wave 
number  which  can  be  adjusted 
for  sensors  of  different  frequen¬ 
cies. 

^r(r)r(2r) 

-  _  2jtpr(2-r)2'-'»'  2- 

r  ( 1  -  r)  T  ( 1  +  r)  “ 

R90  -  pv~  1 
pt)  +  1 

Volume3 

o(a<a  )  -  5^/(1 -g2)Wa 
c  In  (10)  sb 

pucos  a-sb 
pocosa  +  sb 

sb  =  a/i  -o2sin2a 

vol  =  volume  scattering  albedo 

p  =ratio  of  sediment  mass  den¬ 
sity  to  water  mass  density 

v  =  ratio  of  sediment  sound 
speed  to  water  sound  speed 

ccc  =  asin(l/v)  =  surface  critical 
angle 

Composite 

Roughness3 

o  =  4Vjcos4aflV 

W  =  p(2/:asina)  2(r+1) 

Kfr,^n.  )  _  ( (P  —  l)2sin2a  +  p2  +  u-2)2 
(pcosa  +  Jx>~2  -  sin2a) 4 

-  ((P- l)2sin2a  +  p2  +  u-2)2 
((1  -p)2sin2a  +  p2  +  o-2)2 

r  =  y/2-1;  y  =  surface  roughness 
spectrum  exponent 

P  =  surface  roughness  spectral 
strength 

p  =ratio  of  sediment  mass  den¬ 
sity  to  water  mass  density 

v  =  ratio  of  sediment  sound 
speed  to  water  sound  speed 

ac  =  asin(l/v)  =  surface  critical 
angle 

a.  Nayar  et.al.,  1992. 


b.  Torrance  et.  al.,  1967. 


c.  Jackson  et.  al.,  1986. 


Best  Fit  Scattering  Models 


Ridge  Flank 


FIGURE  5.  Best  fit  scattering  models  as  given  in  Table  1  to  the  ridge  flank  data  set. 

tude  grew  proportionally  to  the  range  from  the  sensor  squared  was  added  to  the  synthetically 
generated  images  of  backscatter  and  bathymetry  to  test  the  algorithm  robustness  to  noise.. 


5  Detailed  Local  Algorithm 

The  local  reconstruction  algorithm  reconstructs  the  elevation  profile  of  the  current  ping  in  fish 
coordinates.  This  reconstruction  uses  local  information  to  generate  an  accurate  reconstruction 
and  assumes  that  locally  the  fish  moves  in  a  straight  line 
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5.1  Initial  Estimate  of  Local  Maps 

To  start  the  algorithm,  an  initial  estimate  of  the  elevation  surface  must  be  created  from  the 
sparse  elevation  control  points.  This  initial  estimate  ensures  that  the  starting  point  for  the  sur¬ 
face  is  close  to  its  actual  value.  The  generation  of  the  initial  estimate  in  the  general  case  is 
straightforward.  The  Delaunay  triangulation  of  the  (x,y)  coordinates  of  the  elevation  points  is 
determined  using  well  known  algorithms  from  computational  geometry.  Then  the  elevation 
values  between  control  points  are  determined  by  linearly  interpolating  along  the  plane  fit  to 
each  triangle  in  the  triangulation.  This  interpolation  is  quick  and  sufficient  for  our  algorithm. 
To  simplify  matters  somewhat,  we  have  assumed  that  the  control  points  lie  on  a  grid;  hence, 
computation  of  the  Delaunay  triangulation  is  unnecessary  because  each  comer  in  the  grid 
determines  a  triangle  in  the  triangulation. 

The  cylindrical  nature  of  the  sensor  creates  a  sonar  image  whose  columns  correspond  to 
acoustic  samples  which  can  be  converted  to  range  and  whose  rows  correspond  to  the  position 
of  the  fish  for  each  ping.  Thus  the  image  can  be  parameterized  in  terms  of  the  ping  and  the 
sample  (p,s).  If  the  assumption  is  made  that  locally  the  fish  moves  perpendicularly  to  the 
swath  and  in  a  horizontal  straight  line,  then  locally  the  derivatives  of  the  elevation  surface  can 
be  calculated  by  taking  finite  differences  between  consecutive  rows  and  columns.  The  finite 
differences  along  the  columns  will  give  the  derivative  of  the  surface  with  respect  to  range  zs, 
and  the  finite  difference  taken  in  the  row  direction  will  give  the  derivative  of  the  surface  with 
respect  to  the  distance  travelled  along  the  assumed  straight  line  fish  path  zp.  This  first  interpo¬ 
lation  results  in  three  maps  embedded  in  the  echo  intensity  image  that  correspond  to  the  eleva¬ 
tion  at  each  pixel  z(p,s),  the  derivative  of  the  surface  with  respect  to  the  range  zs(p,s)  and  the 
derivative  of  the  surface  with  respect  to  distance  travelled  zp(p,s). 

Once  the  initial  estimates  for  the  surface  and  its  derivatives  are  generated,  the  local  recon¬ 
struction  for  each  ping  can  begin.  Local  maps  of  elevation  and  its  derivatives  are  created  from 
the  rows  of  the  image  surrounding  the  current  ping.  The  size  of  the  image  (i.e.,  number  of 
rows  included)  can  vary.  Small  images  should  be  used  when  the  position  of  the  fish  changes 
greatly  between  pings;  large  images  should  be  used  when  the  fish  basically  moves  in  a  straight 
line.  Larger  local  images  are  better  because  the  larger  the  image,  the  less  likely  boundary 
effects  will  affect  the  reconstruction.  Once  these  local  maps  are  created,  the  estimation  of  the 
parameters  of  the  scattering  parameters  can  proceed. 

5.2  Local  Scattering  Parameter  Estimation 

The  scattering  model  of  a  surface  is  a  function  of  the  incidence  angle  of  the  signal  striking  the 
surface  and  some  set  of  parameters  p  =  (pl,p2,...)  that  depend  on  the  material  properties  of  the 
surface  and  the  sensor  ct  =  a(a,p).  Even  if  the  incidence  angle  of  the  surface  is  know  at  every 
point  in  the  image,  the  expected  intensity  values  predicted  by  the  scattering  model  cannot  be 
determined  because  the  scattering  parameters  of  the  surface  are  not  known.  However,  given 
the  intensity  image  of  the  seafloor,  the  values  of  the  scattering  parameters  can  be  predicted 
using  least  squares  fitting  of  the  scattering  model  to  the  backscatter  data. 

Because  the  seafloor  is  made  up  of  regions  of  similar  material  type  (which  controls  the  values 
of  the  scattering  model  parameters),  it  is  justifiable  to  assume  that,  in  local  neighborhoods  on 
the  seafloor,  scattering  parameters  will  have  similar  values.  Using  this  assumption,  we  can 
break  the  backscatter  image  into  regions  where  the  scattering  parameter  values  will  be  con- 
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stant.  For  convenience  we  have  made  these  constant  regions  square  and  of  adjustable  size. 
Because  the  scattering  parameters  are  constant  in  square  regions  that  are  larger  than  one  pixel 
they  have  a  lower  resolution  than  the  elevation  map.  Each  point  in  the  image  has  an  equation 
of  the  form  Sj  =  a(aPp),  so  increasing  the  size  of  the  constant  regions  (hence  the  number  of 
points  where  p  is  constant)  will  increase  the  number  of  equations  which  can  be  used  to  solve 
for  the  scattering  parameters  using  non-linear  minimization. 

Under  this  assumption,  in  each  region  the  parameters  p  that  best  fit  the  data  S)  with  error  8j  in 
the  control  point  region  can  be  found  by  minimizing  the  chi  square  merit  function  for  fitting 
the  scattering  model  to  the  data. 


2 

X  (P) 


(EQ  3) 


This  function  was  minimized  with  the  Levenberg-Marquadt  Method  for  non-linear  least 
squares  minimization[27].  The  resulting  parameters  will  be  the  best  fit  parameters  to  the  data 
in  the  region  given  the  scattering  model  and  the  intensity  data. 

A  virtue  of  this  fitting  method  is  that  different  scattering  models  can  be  used  simply  by  chang¬ 
ing  the  function  supplied  to  the  minimization  method.  In  this  way  different  scattering  models 
can  be  tested  on  the  data  to  determine  the  one  that  is  the  most  appropriate.  The  scattering 
model  parameters  are  re-estimated  each  time  a  new  estimate  of  the  incidence  angles  is  deter¬ 
mined  from  the  new  elevation  map.  Convergence  of  the  elevation  map  makes  the  incidence 
angles  converge  to  their  true  values;  as  a  result,  the  estimates  of  the  scattering  parameters 
should  converge  to  their  correct  value.  The  scattering  parameters  are  considered  to  be  close  to 
their  true  values  when  the  value  of  %2  is  close  to  1 .0.  When  this  is  true  for  all  control  point 
regions,  the  scattering  parameters  are  considered  to  be  converged  for  the  surface.  This  is  dem¬ 
onstrated  by  testing  done  on  synthetic  data  where  the  true  albedo  is  known. 

A  shortcoming  of  this  method  for  parameter  estimation  is  that  it  cannot  always  estimate  all  of 
the  scattering  parameters  in  the  scattering  model.  Fortunately,  this  generally  occurs  when  the 
parameter  contributes  little  to  the  form  of  the  scattering  model.  A  specific  example  of  this  dif¬ 
ficulty  occurs  when  trying  to  estimate  o  in  the  Torrance  and  Sparrow  model  at  large  values  of 
a.  This  is  precisely  because,  at  large  values  of  a,  the  specular  component  of  scattering  is  neg¬ 
ligible,  so  numerically  a  correct  estimate  of  a  is  not  easy  to  obtain. 

It  is  also  difficult  to  estimate  parameters  correctly  when  there  is  a  strong  coupling  between 
parameters.  This  coupling  means  there  exist  multiple  sets  of  parameters  that  generate  the 
same  scattering  model  form.  The  data  will  fit  the  scattering  model  equally  well  for  multiple 
sets  of  parameters,  so  the  minimization  method  has  no  way  to  select  the  correct  scattering 
parameters.  However,  the  form  of  the  scattering  model  will  generally  be  estimated  correctly. 
An  example  of  this  coupling  occurs  in  the  diffuse  model  where  p  and  £  are  multiplicatively 
coupled  in  the  roughness  term  of  the  model. 

Another  difficulty  in  estimating  the  scattering  parameters  is  that  some  parameters  have  defi¬ 
nite  constraints  on  the  values  that  they  can  take.  For  instance  all  of  the  parameters  cannot 
assume  negative  values.  Occasionally,  if  the  data  is  noisy  or  the  surface  estimate  is  very 
wrong,  the  estimate  of  the  parameters  creates  a  form  for  the  scattering  model  that  fits  the  data 
well  and  yet  is  generated  by  parameters  that  are  not  physically  feasible.  To  circumvent  this 
problem,  we  add  terms  to  the  scattering  model  in  the  chi  square  fitting.  These  terms  become 
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(B) 
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(C)  Initial  Backscatter  vs.  Incidence  Angle 


(D)  Final  Backscatter  vs.  Incidence  Angle 


FIGURE  6.  Convergence  of  the  albedo  parameter  for  a  synthetic  data  set  following  the  lambertian 
scattering  model.  (A)  Initial  albedo  isoplot.  (B)  Final  albedo  isoplot.  (C)  Initial  graph  of 
backscatter  vs.  incidence  angle.  (D)  Final  graph  of  backscatter  vs.  incidence  angle. 


very  large  when  the  parameters  go  outside  of  their  expected  ranges.  Thus  when  a  parameter  is 
outside  its  range,  the  value  of  %2  becomes  very  large.  This  will  force  the  parameters  to  be 
adjusted  in  such  a  way  as  to  decrease  %2  and  will  occur  only  when  the  parameter  that  was  out¬ 
side  of  its  range  comes  back  into  range. 

- - - — - — - - - — - - - 
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Each  time  that  the  surface  is  reestimated,  the  scattering  parameters  are  reestimated  as  well.  In 
this  way,  the  scattering  parameters  converge  to  their  best  fit  form  as  long  as  the  surface  is  con¬ 
verging  to  its  true  value.  Figure  6  shows  the  convergence  of  the  albedo  parameter  for  a  syn¬ 
thetic  data  generated  from  a  cylinder.  Initially  the  estimates  of  albedo  for  the  quite  different 
from  there  true  value  of  0.5.  However,  as  the  algorithm  proceeds,  the  estimates  of  the  albedo 
improve  drastically.  The  isoplots  show  that  the  estimates  of  albedo  have  a  much  larger  spread 
and  are  far  from  the  actual  value  of  albedo  when  the  algorithm  starts,  but  at  the  end  of  the 
algorithm,  the  estimates  have  a  small  spread  and  are  close  to  the  expected  value  for  the 
albedo.  The  graphs  show  that  initially  the  estimates  of  incidence  angle  are  incorrect,  but  at  the 
end  of  the  iterations  to  improve  the  shape  of  the  surface,  the  incidence  angle  are  estimated 
correctly,  so  the  albedo  estimates  are  much  more  accurate.  This  test  on  synthetic  data  demon¬ 
strates  that  the  algorithm  for  estimating  scattering  parameters  is  feasible. 


5.3  Determining  the  Elevation  Map 

The  backscatter  image  of  the  seafloor  does  not  contain  enough  information  about  the  shape  of 
the  seafloor  to  determine  it  exactly.  If  we  had  exact  knowledge  of  the  scattering  model  for  the 
surface  as  well  as  its  parameter,  we  would  still  find  it  impossible  to  determine  the  shape  of  the 
seafloor  uniquely.  Each  pixel  in  the  backscatter  image  can  be  directly  related  to  the  angle  of 
incidence  of  the  acoustic  pulse  with  the  surface,  but  the  angle  of  incidence  is  a  function  of  the 
elevation  of  the  surface  and  its  two  surface  derivatives.  Thus  with  perfect  scattering  informa¬ 
tion  for  every  pixel  there  will  be  only  one  equation  for  every  3  unknowns  S  =  a(a(z,zs,zp),p) 
In  order  to  solve  for  the  elevation  map,  some  assumptions  must  be  made. 

Our  approach  is  to  solve  for  the  surface  by  minimizing  some  global  criteria  that  measure  the 
goodness  of  fit  to  the  data.  Because  the  data  can  fit  an  infinite  number  of  surfaces,  constraints 
must  be  imposed  on  the  shape  of  the  surface  to  decrease  the  total  number  of  possible  surfaces 
to  one.  The  procedure  is  to  define  an  error  functional  that  strikes  a  balance  between  the  good¬ 
ness  of  fit  of  the  surface  to  the  data  and  the  need  to  constrain  the  surface  so  that  a  unique  solu¬ 
tion  can  be  obtained.  Minimization  of  such  an  error  functional  using  shape  from  shading 
techniques  appears  frequently  in  computer  vision.  After  trying  many  different  forms  for  the 
error  functional,  we  determined  that  the  form  given  in  Equation  4  gave  the  best  result  for  the 
maps  of  z,  zs  and  zp.  This  error  functional  is  similar  to  that  developed  by  Horn  and  Brooks  for 
shape  from  shading  with  an  integrability  penalty  term  [15]  except  that  the  seafloor  backscatter 
model  is  a  function  of  elevation  as  well  as  its  derivatives. 


\\(S(s,p)  -a{zfz^))2dsdp*]\^zt-^  +(zp-fr)  dsdp  +  l.\\[ 


2  2  2  2 
+  Z sp  +  *"ps  +  Zpp 


dsdp 


(EQ4) 


The  first  integral  in  this  functional  ensures  that  the  data  matches  the  scattering  model  in  the 
least  squares  sense  and  the  second  integral  is  a  constraint  term  that  tries  to  maintain  the  inte¬ 
grability  of  the  surface.  These  two  integrals  are  sufficient  to  uniquely  determine  the  surface 
given  the  intensity  scattering  data;  however,  the  solution  obtained  by  using  only  these  two 
terms  will  diverge  if  too  much  noise  is  present  in  the  image.  The  third  integral  ensures  that  the 
solution  does  not  diverge  by  requiring  that  the  reconstructed  elevation  map  vary  smoothly. 
The  smoothness  term  will  force  the  estimates  of  the  surface  to  remain  close  together,  thus  pre¬ 
venting  transient  errors  in  the  estimation  from  causing  the  surface  estimates  to  diverge.  The 
smoothness  of  the  seafloor  will  vary  from  one  geological  province  to  the  next,  so  the  smooth- 
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ness  of  the  reconstructed  surface  can  be  adjusted  by  changing  the  weight  X  on  the  smoothness 
error  term.  Similarly,  the  degree  to  which  integrability  is  enforced  on  the  surface  is  influenced 
by  the  value  of  ji  in  the  functional.  These  values  were  empirically  determined  to  give  the  best 
results  when  X  =  0.1  and  |A  =  0.01. 

The  error  functional  is  minimized  by  applying  the  calculus  of  variations  to  determine  the 
Euler  equations  for  the  functional.  Inserting  finite  differences  estimates  for  the  continuous 
derivatives  will  discretize  the  Euler  equations.  These  three  discrete  Euler  equations  can  then 
be  solved  for  the  three  surface  unknowns  z,  zs  and  zp  to  get  the  three  iterative  update  equations 
given  in  Equation  5. 


n+\  |i  z"(i,j+  1)  -zn(i,j)  1  ..  (  n  n  njjda  X  n  , . 
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The  bar  represents  an  average  of  the  four  nearest  neighbors  of  a  point.  The  scattering  model  is 
a  function  of  the  zs,  zp  and  z  because,  as  is  shown  in  the  appendix,  the  incidence  angle  can  be 
defined  as  given  in  Equation  6. 
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5.4  Algorithm 

By  iteratively  adjusting  the  elevation  map,  the  elevation  derivative  maps,  and  the  scattering 
parameter  maps  until  the  estimates  change  very  little,  the  algorithm  will  find  the  best  estimate 
of  the  elevation  in  fish  coordinates  given  the  data.  The  complete  local  surface  reconstruction 
algorithm  is  as  follows. 

For  all  pings  in  the  acoustic  backscatter  image: 

1 .  Linearly  interpolate  z  between  the  known  elevation  control  points. 

2.  Assuming  that  locally  the  fish  moves  in  a  straight  line,  create  local  derivative  maps  cen¬ 
tered  around  the  current  ping  by  taking  finite  differences  between  rows  and  columns  in 
the  interpolated  elevation  image. 

3.  Estimate  the  scattering  parameters  by  fitting  the  data  from  square  regions  in  the  local 
map  to  the  scattering  model  using  non-linear  minimization. 

4.  Determine  the  new  values  for  zs,  zp  and  z  using  Equation  5. 

5.  Repeat  steps  3  and  4  until  the  intensity  predicted  by  the  reconstructed  surface  is  close  to 
the  intensity  given  by  the  data  and  the  elevation  estimate  do  not  change  very  much. 
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5.5  Local  Results 

The  convergence  of  the  elevation  estimate  for  a  synthetically  generated  cylinder  that 
employed  the  lambertian  scattering  model  and  had  elevation  control  points  spaced  every  40 
pixels  can  be  seen  in  Figure  7.  A  cross-section  of  the  three  isoplots  at  row  25  shows  the  con¬ 
vergence  of  the  estimate  of  the  elevation  z.  For  this  experiment  on  synthetic  data  a  fixed 
boundary  condition  is  enforced  to  prevent  the  surface  from  diverging  too  greatly  from  its  ini¬ 
tial  estimate.  The  fixed  boundary  will  not  affect  the  final  reconstruction  because  for  local 
reconstruction  we  are  only  concerned  with  the  middle  of  the  surface.  Initially  the  surface  is 
flat,  but  at  the  end  of  the  iterations  the  surface  takes  on  the  shape  of  the  curved  cylinder.  This 
test  on  synthetic  data  shows  that  the  surface  refinement  algorithm  can  use  backscatter  data  to 
improve  the  estimate  the  surface  of  the  seafloor  given  by  sparse  elevation  data  points.  The 
convergence  of  the  local  elevation  estimate  on  a  set  of  real  data  is  shown  in  Figure  8.  The  ini¬ 
tial  estimate  for  the  real  surface  was  generated  from  points  lying  at  the  comers  of  a  250  m  by 
150  m  rectangle.  Each  pixel  in  the  images  represents  an  area  of  0.5  m  by  0.5  m.  The  sequence 
of  isoplots  shows  the  initial  flat  estimate,  the  final  reconstructed  surface  and  the  actual  surface 
predicted  from  the  bathymetry  given  by  the  sensor  for  a  small  section  of  the  reconstructed  sur¬ 
face.  The  real  surface  is  not  recovered  exactly,  but  the  final  estimate  is  much  better  than  the 
initial  estimate.  The  graph  of  a  row  from  each  of  the  three  isoplots  shows  that  the  final  esti¬ 
mate  of  the  surface  is  a  much  better  estimate  of  the  surface  that  the  initial  estimate.  Figure  9 
shows  the  error  histograms  for  difference  between  the  true  surface  and  the  initial  estimate  and 
the  true  surface  and  the  final  estimate.  The  histograms  clearly  show  that  the  error  decreases 
from  initial  estimate  of  the  surface  to  final  estimate  of  the  surface.  Figure  9  also  shows  the 
average  percentage  decrease  in  the  error  between  the  actual  surface  and  the  estimates  surface. 
The  error  decreases  monotonically  to  a  local  minimum.  The  final  elevation  estimate  does  nor 
completely  agree  with  the  actual  surface  because  the  scattering  data  does  not  correspond 
exactly  to  that  predicted  by  the  surface  and  the  surface  generate  must  remain  smooth,  which 
will  prevent  any  rapid  changes  in  the  curvature  of  the  surface. 


6  Local  to  Global 

Acoustic  images  can  be  deceiving  because  they  are  created  by  stacking  consecutive  pings  on 
top  of  each  other.  This  stacking  has  the  effect  of  tricking  the  viewer  into  thinking  that  the 
image  was  created  by  towing  the  fish  in  a  horizontal  straight  line.  In  general,  the  fish’s  path 
will  deviate  from  a  straight  line  as  it  is  towed  through  the  water  due  to  variable  tension  in  the 
tether,  surface  turbulence  and  water  currents.  Since  the  pose  of  the  fish  influences  the  back¬ 
scatter  that  the  sensor  detects  (because  the  angle  of  incidence  of  the  acoustic  pulse  is  a  func¬ 
tion  of  the  surface  normal  and  the  position  of  the  fish)  the  backscatter  will  have  meaning  only 
in  the  coordinates  of  the  fish  for  the  current  ping.  In  the  previous  section,  a  method  for  recon¬ 
structing  the  elevation  profile  for  the  current  ping  in  current  fish  coordinates  was  presented. 
However,  these  local  reconstructions  are  useless  unless  they  are  transformed  to  some  common 
coordinate  system  where  true  global  proximity  can  be  determined.  Given  this  proximity, 
meaningful  2  1/2  D  elevation  maps  of  the  seafloor  can  be  generated.  Transformation  of  the 
elevation  profiles  for  each  ping  to  global  coordinates  is  straightforward,  while  transformation 
of  acoustic  backscatter  which  is  necessary  for  global  elevation  refinement  is  not. 
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FIGURE  7.  Convergence  of  local  elevation  for  a  synthetic  image  of  a  cylinder.  (A)  Initial  Estimate. 

(B)  Final  reconstruction.  (C)  Actual  surface.  (D)  Row  25  from  each  preceding  isoplot. 

6.1  Binning  and  the  3D  Transformation 

The  global  coordinates  are  arbitrarily  chosen  to  be  the  coordinates  of  the  fish  for  the  first  ping 
in  the  image.  Given  that  the  fish  moves  in  a  path  that  is  close  to  a  horizontal  straight  line,  this 
assignment  will  make  the  transformation  from  local  to  global  coordinates  small,  hence  mak¬ 
ing  the  global  surface  appear  similar  to  the  local  reconstructions.  The  global  elevation  surface 
will  be  generated  on  a  square  grid  whose  origin  corresponds  to  the  origin  of  the  global  coordi¬ 
nate  system.  The  rows  of  the  elevation  map  correspond  to  the  y  global  coordinate  and  the  col¬ 
umns  correspond  to  the  x  coordinate. 
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Elevation  for  Ridge  Flank 


(row  55) 


FIGURE  8.  Convergence  of  the  local  elevation  estimate  for  a  small  area  of  ridge  flank.  (A)  Initial 
estimate.  (B)  Final  reconstructed  surface.  (C)  Actual  surface.  (D)  Row  55  from  each  of 
the  three  preceding  isoplots. 


For  each  ping,  the  local  elevation  profiles  which  are  in  cylindrical  fish  coordinates  are  trans¬ 
formed  into  cartesian  global  coordinates  using  the  pose  of  the  fish  as  is  shown  in  Equation  7. 
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Local  Elevation  Error  Histograms  ,D.  Average  Elevation  Error  Convergence 

(A)  (B) 

ridge  flank  ridge  flank 


FIGURE  9.  (A)  Initial  and  final  elevation  error  histograms  for  ridge  flank  and  (B)  percent  decrease 

in  average  error  as  a  function  of  algorithm  iterations. 


The  x  and  y  global  coordinates  of  a  point  determine  the  grid  cell  in  the  global  map  in  which  to 
place  the  point.  Depending  on  the  resolution  of  the  global  map,  multiple  points  may  fall  into 
the  same  grid  cell.  If  this  is  the  case,  the  global  z  values  of  all  of  the  points  are  averaged  to  get 
the  initial  estimate  of  the  global  elevation  for  that  grid  cell.  In  general,  the  resolution  of  the 
global  map  is  set  such  that,  on  average,  each  grid  cell  will  have  one  local  point  fall  into  it.  This 
resolution  will  strike  a  balance  between  the  gaps  and  overlap  that  might  appear  in  the  global 
map.  The  transformation  and  binning  of  the  local  elevation  estimates  creates  the  initial  global 
surface  estimate  for  the  that  will  be  refined  again  using  the  backscatter  data.  The  problem  then 
is  to  find  a  way  to  transform  the  backscatter  data  such  that  the  data  is  still  useful  and  little  or 
none  of  its  information  is  lost. 

Transforming  and  then  averaging  of  the  local  elevation  points  that  fall  in  the  same  grid  cell  is 
appropriate  because  the  addition  of  elevation  values  is  meaningful.  However,  a  meaningful 
way  to  combine  backscatter  is  not  apparent.  The  backscatter  data  is  a  function  of  the  position 
of  the  sensor  and  can  vary  over  the  points  that  fall  in  the  same  grid  cell.  It  is  entirely  possible 
for  scattering  values  that  fall  in  the  same  grid  cell  to  originate  from  pings  whose  fish  positions 
are  very  different.  Simple  averaging  of  these  types  of  scattering  values  will  be  useless  because 
the  position  of  the  fish  of  the  averaged  scattering  value  is  undefined.  The  proposed  solution  to 
this  problem  is  to  create  a  bin  data  structure  that  stores  all  of  the  backscatter  values  and  the 
corresponding  position  of  the  fish  when  the  data  was  collected  for  points  that  fall  in  each  grid 
cell.  This  solution  ensures  that  none  of  the  backscatter  data  is  lost  and  that  the  backscatter  data 
can  be  used  because  the  position  of  the  fish  is  known,  through  the  ping  that  falls  in  the  grid 
cell. 


6.2  Importance  of  Local  Reconstructions 

The  local  reconstructions  were  necessary  before  the  transformation  to  global  coordinates 
because  of  the  cylindrical  nature  of  the  data  in  fish  coordinates.  The  estimate  of  local  eleva- 
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tion  determines  the  local  x  and  z  coordinates  of  the  data  point,  so  if  the  elevation  estimate  is 
bad,  these  coordinates  will  be  bad;  the  data  point  could  be  transformed  into  the  wrong  global 
grid  cell.  Since  the  grid  cell  into  which  the  point  is  transformed  remains  fixed  after  transfor¬ 
mation,  this  error  could  mar  the  global  reconstruction.  The  better  the  local  elevation  estimate 
is  the  better  the  global  reconstruction  will  be.  If  the  local  reconstructions  are  perfect,  then  no 
global  reconstruction  needs  to  be  done  once  the  points  are  transformed. 


7  Detailed  Global  Algorithm 

The  local  estimates  are  transformed  into  global  coordinates  and  averaged  in  bins  to  create  the 
initial  estimate  for  the  global  elevation  surface.  The  estimates  of  the  derivatives  of  the  surface 
can  then  be  derived  by  taking  the  finite  differences  of  the  elevation  surface  along  the  rows  and 
columns  of  the  map.  These  estimates  will  be  refined  into  the  final  estimate  for  the  shape  of  the 
seafloor  in  global  coordinates  given  the  backscatter  data.  These  estimates  will  also  be  used  to 
determine  the  final  scattering  parameter  maps.  The  global  reconstruction  algorithm  is  quite 
similar  to  the  local  one,  but  takes  into  account  the  cartesian  form  of  the  elevation  map  and  the 
binning  of  the  scattering  data. 


7.1  Global  Scattering  Parameter  Estimation 

Once  the  initial  surface  estimates  are  generated,  estimates  for  the  scattering  parameters  are 
generated.  Just  as  in  the  local  algorithm,  the  global  map  is  broken  into  square  regions  where 
the  scattering  parameters  are  considered  to  be  fixed.  This  assumption  then  allows  solution  for 
the  scattering  parameters  given  the  data  using  non-linear  minimization.  The  difference  from 
the  local  algorithm  is  that  the  data  in  the  regions  might  come  from  different  places  in  the  orig¬ 
inal  image  and  that  each  grid  cell  could  contain  more  than  one  backscatter  data  point.  How¬ 
ever,  the  position  of  the  fish  that  produced  each  data  point  is  stored  in  the  bin  data  structure,  so 
we  can  use  this  information  when  doing  the  minimization.  In  effect,  all  that  has  changed  from 
the  local  algorithm  is  the  calculation  of  the  incidence  angle  for  each  data  point  and  the  fact 
that  we  must  use  all  of  the  data  points  that  fall  in  each  grid  cell.  To  get  the  best  fit  scattering 
parameters  given  the  backscatter  data,  we  minimize  the  following  merit  function. 
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In  Equation  8  all  of  the  data  is  considered  because  all  of  the  bins  in  the  region  and  all  of  the 
data  in  the  bins  are  used  in  the  summation.  The  surface  normal  n  is  fixed  for  each  bin  while 
the  position  of  the  fish  P  varies  over  the  bins  and  the  data  in  the  bins,  so  the  incidence  angle  a 
is  predicted  correctlygiven  the  surface  estimate  for  each  data  point.  The  functional  form  of  the 
scattering  model  does  not  change  between  local  and  global  algorithms,  so  this  merit  function 
can  be  minimized  using  the  Levenberg-Marquadt  method  as  was  done  in  the  local  algorithm. 


7.2  Global  Elevation  Reconstruction 

Given  the  estimates  of  the  scattering  parameters  in  the  global  maps,  the  next  stage  in  the  glo¬ 
bal  algorithm  is  to  adjust  the  estimates  of  the  seafloor  so  that  they  fit  the  data  according  to 
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some  global  error  functional.  As  was  done  for  the  local  algorithm,  we  have  constructed  an 
error  functional  that  measures  the  goodness  of  fit  of  the  data  to  the  surface  estimate  while 
maintaining  the  integrability  of  the  surface  and  the  convergence  of  the  algorithm.  The  func¬ 
tional  given  in  Equation  9  is  similar  to  the  local  error  functional  except  that  all  of  the  data  in 
each  bin  is  considered.  This  functional  is  similar  to  those  used  in  shape  from  shading  tech¬ 
niques  where  images  of  the  scene  from  multiple  cameras  are  available  [14]. 
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The  first  integral  compares  all  of  the  data  to  the  scattering  model  estimate  which  is  a  function 
of  the  current  surface  estimates  as  well  as  the  position  of  the  fish  when  the  data  was  recorded. 
The  second  integral  maintains  integrability  of  the  surface  and  the  third  integral  keeps  the  sur¬ 
face  estimates  from  diverging.  Iterative  update  equations  are  generated  using  the  calculus  of 
variations  and  discrete  forms  for  the  derivatives  in  the  exact  same  way  as  was  done  for  the 
local  algorithm. 

The  control  of  the  global  algorithm  is  identical  to  the  local  algorithm.  Initial  estimates  are 
generated;  then  the  estimates  for  scattering  parameters,  the  surface  and  its  derivatives  are  iter¬ 
atively  refined  until  the  surfaces  converge.  These  estimates  are  the  best  global  estimates  of  the 
elevation  map  of  the  seafloor  and  its  corresponding  parameter  maps  given  the  acoustic  back- 
scatter  image.  The  global  algorithm  for  refinement  of  the  seafloor  “averages”  the  scattering 
data  in  each  bin  during  the  reconstruction,  so  no  information  is  lost. 

8  Global  Results 


The  global  reconstruction  process  was  first  tested  on  synthetic  data  to  ensure  its  viability. 
Figure  10  shows  various  stages  in  the  global  reconstruction  process  on  a  set  of  synthetic  data 
generated  as  the  fish  moved  in  a  semi-circular  path  around  a  cylinder.  The  initial  estimate  for 
the  surface  in  no  way  resembles  a  cylinder  because  the  shape  of  the  cylinder  is  distorted  by  the 
curved  path  of  the  fish.  The  final  local  elevation  map  shows  a  surface  with  greater  curvature, 
but  the  geometry  of  the  surface  is  still  extremely  skewed.  Not  until  the  surface  is  transformed 
to  global  coordinates  in  the  initial  global  elevation  map,  does  the  shape  of  the  cylinder 
become  apparent.  The  curvature  generated  by  the  local  reconstructions  is  present  in  the  initial 
global  elevation  map,  but  the  true  curvature  of  the  cylinder  is  not  realized  until  global  recon¬ 
struction  algorithm  completes.  The  final  global  reconstruction  has  the  same  shape  as  the 
actual  cylinder  imaged  except  at  the  boundaries  where  the  elevation  was  kept  fixed  to  con¬ 
strain  the  surfaces  generated.  As  is  shown  in  the  density  of  overlap  image,  the  path  of  the  fish 
for  this  synthetic  data  set  caused  the  highest  density  (white)  of  overlapping  data  to  occur  at  the 
center  of  the  fish’s  circular  path.  This  synthetic  data  set  clearly  shows  how  a  nonlinear  fish 
path  can  change  the  geometry  of  the  surface  imaged,  and  it  showed  that  our  global  algorithm 
works  to  reconstruct  the  actual  geometry  of  the  surface  in  a  global  sense  given  the  path  of  the 
fish. 


22 


FIGURE  10.  Convergence  of  the  global  elevation  estimate  for  a  synthetic  data  set  where  the  fish  moves 
in  a  circle  about  a  cylindrical  surface.  (A)  Initial  local  estimate.  (B)  Final  local  estimate. 
(C)  Initial  global  estimate.  (D)  Final  global  estimate.  (E)  Actual  surface  and  path  of  sonar 
fish.  (F)  Overlap  of  data. 


The  result  of  the  global  reconstruction  process  performed  on  the  ridge  flank  data  set  using  the 
lambertian  scattering  model  is  given  in  Figure  11.  The  initial  estimate  of  the  surface  was  gen¬ 
erated  from  a  regular  grid  of  elevation  control  points  spaced  every  40  meters.  The  recon¬ 
structed  corresponds  very  well  to  the  actual  surface  although  some  sharp  features  are  lost  in 
the  reconstruction  due  to  the  smoothing  term  in  the  elevation  error  functional.  The  recovered 
albedo  shows  that  the  albedo  increases  perpendicularly  to  the  direction  of  travel  of  the  fish  and 
that  in  general  the  albedo  of  the  surface  is  very  low  which  is  in  agreement  with  the  soft  mate¬ 
rials  that  made  up  the  ridge  flank.  Verification  of  the  recovered  albedo  is  not  possible  because 
there  are  no  absolute  measurements  of  scattering  parameters.  The  overlap  image  (white  corre¬ 
sponds  to  a  high  degree  of  overlap)  shows  that  the  fish  moved  in  a  fairly  straight  line,  how¬ 
ever,  there  do  exist  places  in  the  image  where  scan  lines  overlapped. 
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FIGURE  11.  Global  reconstruction  on  the  ridge  flank  data  set.  (A)  Recovered  elevation  map.  (B) 

Actual  global  elevation  surface  given  path  of  fish.  (C)  Recovered  albedo  parameter.  (D) 
Overlap  of  data  in  global  image. 
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The  estimation  of  the  scattering  model  parameters  from  the  reconstructed  elevation  maps  can 
be  used  to  segment  the  seafloor  based  on  scattering  model  parameters.  Because  scattering 
model  parameters  depend  on  the  material  properties  of  the  seafloor,  this  segmentation  essen¬ 
tially  finds  regions  of  the  seafloor  of  similar  material.  Figure  12  depicts  how  the  seafloor  can 
be  segmented  based  on  material  type.  The  data  was  generated  by  imaging  an  axial  valley 
region  of  the  seafloor.  The  reconstruction  algorithm  determined  that  there  existed  regions  of 
high  and  low  albedo  in  the  image  when  the  lambertian  scattering  model  was  applied.  In  the 
axial  valley,  these  regions  could  come  from  the  sandy  bottom  of  the  valley  close  to  the  fish 
(with  a  low  albedo)  and  the  rocky  sides  of  the  valley  far  from  the  fish.  This  result  shows  that 
qualitative  segmentation  of  the  seafloor  is  possible  with  our  reconstruction  and  parameter  esti¬ 
mation  algorithm. 

Estimation  of  the  scattering  model  parameters  using  the  volume  scattering  model  was  con¬ 
ducted  on  data  set  taken  of  the  ridge  flank  at  a  different  location  from  the  one  used  in  the  pre¬ 
vious  figures.  The  three  isoplots  of  the  estimated  parameters  along  with  the  corresponding 
actual  surface  and  bathymetric  data  are  given  in  Figure  13.  The  volume  scattering  albedo  is 
generally  small  which  is  in  agreement  with  the  previous  ridge  flank  data  of  Figure  4.  The  two 
density  ratios  are  within  the  expected  ratios  (greater  than  one  and  generally  less  than  three). 


FIGURE  12.  Segmentation  of  data  set  by  albedo  parameter.  (A)  Actual  elevation  of  surface.  (B) 

Estimate  albedo  isoplot.  (C)  Estimated  albedo  image  showing  regions  of  high  and  low 
albedo  giving  a  segmentation  of  the  seafloor  based  on  material  type. 


25 


9  Conclusion  and  Contributions 


In  this  paper  we  have  presented  an  algorithm  that  generates  global  maps  of  elevation  and  scat¬ 
tering  model  parameters  from  a  side  scan  sonar  backscatter  image,  sparse  bathymetric  data  of 
the  seafloor  and  the  path  of  the  sensor  through  the  water.  This  algorithm  works  independently 
of  the  scattering  model  used  and,  in  regions  bounded  by  known  elevation  values,  determines 
the  best  set  of  scattering  parameters  given  the  scattering  model  using  non-linear  least  squares 
fitting.  The  elevation  map  of  the  seafloor  is  determined  by  minimizing  a  global  error  func¬ 
tional  that  measures  the  error  between  the  scattering  data  and  the  backscatter  predicted  by  the 
elevation  map  while  taking  into  account  the  deviations  of  the  path  of  the  fish  from  a  straight 
line.This  algorithm  required  the  development  of  a  sophisticated  method  for  merging  backscat¬ 
ter  data  in  a  global  coordinate  system  that  minimizes  the  loss  of  information. 

The  algorithm  has  been  shown  to  correctly  predict  elevation  and  scattering  parameter  values 
for  synthetic  and  real  data  using  different  scattering  models.  This  flexibility  allows  the  algo¬ 
rithm  to  be  used  by  many  different  types  of  sidescan  sensors  and  on  seafloors  of  differing 
material  properties.The  algorithm  has  also  been  shown  to  build  accurate  elevation  maps  and 
segment  scattering  parameters  where  appropriate  on  real  data 

The  many  stages  of  the  algorithm,  including  the  nonlinear  minimization  in  regions  of  constant 
scattering  parameters  and  the  global  relaxation,  readily  lend  themselves  to  a  parallel  approach 
for  the  algorithm.  The  speed  of  the  algorithm  could  also  be  increased  by  implementing  a 
coarse  to  fine  version  of  this  algorithm.  Given  the  increase  in  run  time  that  a  parallel  and 
coarse  to  fine  approach  would  provide,  it  is  feasible  that  this  algorithm  could  eventually  be 
implemented  as  a  real-time  mapping  and  navigational  package  for  an  autonomous  underwater 
vehicle. 
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12  Appendix:  Local  Angle  of  Incidence 


The  surface  F  that  is  going  to  be  reconstructed  can  be  parameterized  in  R  and  y  in  the  follow¬ 
ing  coordinates. 


F(R,y) 


x(R,y) 
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(EQ  10) 


Assuming  that  the  surface  is  explicit,  at  any  point  (R,y)  the  normal  to  the  surface  is  given  by 
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Differentiating  Equation  10  with  respect  to  R  and  y  gives 
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where  zR  and  zy  are  the  derivatives  of  the  surface  z  in  the  R  and  y  directions  respectively,  so 
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Figure  2  shows  that  for  a  sidescan  sonar 

R  n  =  ||/f|[  cosoc  . 
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So  Figure  14  can  be  solved  for  a  to  get 
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